### Replication code for:
### "Driving While Unauthorized: Auto Insurance in Remains Unchanged when Providing Driver Licenses to Unauthorized Immigrants in California"
### Hans Lueders and Micah Mumper
### Journal of Risk and Insurance
### April 10, 2022




# Introductory remarks ----------------------------------------------------

### This file replicates all figures and tables reported in the main text and appendix of the paper.
### For questions, please reach out to hlueders@stanford.edu or hans.lueders89@gmail.com.

### The replication data contains the following files:
# - d.county: county-year panel data with information on exposure to AB60, all outcomes, and control variables
# - d.aps: data from the Auto Premium Survey. Two observations for each county-year: premiums for experienced and inexperienced drivers.
# - separate files with the results of the power analysis. See separate STATA do-file on how these values were obtained.

### Additionally, replication requires a shapefile of California's counties, which is provided separately.

### The d.county dataset contains the following variables:
# - county: county name
# - year: year
# - yearEVENT: year relative to last pre-treatment year (2014); used for event study plots
# - uninsured: estimated share of uninsured vehicles
# - ClaimsNormStandUMBIBasic: Uninsured motorist BI claims (per 1,000 vehicles)
# - ClaimsNormStandUMBIBasic_log: Uninsured motorist BI claims (per 1,000 vehicles), natural logarithm
# - ClaimsNormNStandBIBasic: Claims volume
# - ClaimsNormNStandBIBasic_log: Claims volume, natural logarithm
# - AvgPremiumStandBIBasic: Average monthly premium for bodily injury
# - AvgPremiumStandBIBasic_log: Average monthly premium for bodily injury, natural logarithm
# - premium1both: Annual premium for inexperienced drivers, both genders
# - premium1both_log: Annual premium for inexperienced drivers, both genders, natural logarithm
# - premium7both: Annual premium for experienced drivers, both genders, natural logarithm 
# - premium7both_log: Annual premium for experienced drivers, both genders, natural logarithm 
# - autos: Number of registered cars
# - autos_log: Number of registered cars, natural logarithm
# - PremiumStandBIBasic: Premium volume
# - PremiumStandBIBasic_log: Premium volume, natural logarithm
# - ExposuresStandBIBasic: Policy volume
# - ExposuresStandBIBasic_log: Policy volume, natural logarithm
# - claimfrequency: Claim frequency (claim volume divided by policy volume)
# - UMBI_BI_ratio: UMBI/BI ratio
# - ClaimsNormStandUMPDBasic: Uninsured motorist PD claims (per 1,000 vehicles)
# - AvgPremiumStandPDBasic: Average monthly premium for property damage
# - AvgPremiumStandPDBasic_log: Average monthly premium for property damage, natural logarithm
# - treat: dummy for post-treatment period
# - exposure: estimated exposure to AB60: share of AB60 licenses issued in 2015
# - exposure_mean: dichotomous indicator of all counties above mean exposure
# - exposure_median: dichotomous indicator of all counties above median exposure
# - exposure_tercile: trichotomous indicator of whether a county is in the first, second, or third tercile on the distribution of exposure
# - exposure_tercileDICH: dichotomous indicator of whether a county falls in the first or second (0) or third tercile on the distribution of exposure
# - exposureNEW: alternative estimated exposure to AB6-: share of AB60 licenses issued in 2015 and 2016
# - exposureNEW_tercile: trichotomous indicator of whether a county is in the first, second, or third tercile on the distribution of exposure (measured in 2015/16)
# - undoc_share: estimated share of county population that is unauthorized
# - undoc_share_tercile: trichotomous indicator of whether a county is in the first, second, or tercile on the distribution of share unauthorized
# - population: county population
# - pop.mean: average population over the sample period
# - popdens_log: log population density
# - income_log: log household income
# - unemprate: unemployment rate

### The d.aps dataset contains the following variables:
# - county: county name
# - year: year
# - yearEVENT: year relative to last pre-treatment year (2014); used for event study plots
# - premium: average premium
# - premium_log: average premium, natural log
# - treat: dummy for post-treatment period
# - treat2: dummy for profile for inexperienced drivers
# - exposure: estimated exposure to AB60: share of AB60 licenses issued in 2015
# - exposure_mean: dichotomous indicator of all counties above mean exposure
# - exposure_median: dichotomous indicator of all counties above median exposure
# - exposure_tercile: trichotomous indicator of whether a county is in the first, second, or third tercile on the distribution of exposure
# - exposure_tercileDICH: dichotomous indicator of whether a county falls in the first or second (0) or third tercile on the distribution of exposure
# - exposureNEW: alternative estimated exposure to AB6-: share of AB60 licenses issued in 2015 and 2016
# - exposureNEW_tercile: trichotomous indicator of whether a county is in the first, second, or third tercile on the distribution of exposure (measured in 2015/16)
# - undoc_share: estimated share of county population that is unauthorized
# - undoc_share_tercile: trichotomous indicator of whether a county is in the first, second, or tercile on the distribution of share unauthorized
# - pop.mean: average population over the sample period
# - popdens_log: log population density
# - income_log: log household income
# - unemprate: unemployment rate

### The files with the results of the power analysis contain the following variables:
# - effect_size: pseudo-treatment effect
# - significant: share of draws that were significant (i.e., power estimate)
# - b: average estimated effect size















# Preparation -------------------------------------------------------------

### Load packages
library(foreign)
library(readxl)
library(lfe)
library(sandwich)
library(lmtest)
library(survminer)
library(ggplot2)
library(gridExtra)
library(xtable)
library(stargazer)
library(survival)
library(lubridate)
library(MASS)

### Set working directory
setwd("") #include file path here

### clear environment
rm(list=ls())

### load data
load("LuedersMumperJRI2022_replicationdata.RData")







# Main text ---------------------------------------------------------------


### Figure 1: Estimated exposure to AB60
# Figure 1a: Exposure by county
shp <- readOGR("shapefiles/CA_counties.shp")
shp.fort <- fortify(shp, region="NAME")
shp.fort <- merge(shp.fort, d.county[d.county$year==2015, c("county", "exposure")], by.x="id", by.y="county", all.x=T)
shp.fort <- shp.fort[order(shp.fort$order),]
shp.fort$exposure_cat <- NA
shp.fort$exposure_cat[shp.fort$exposure >= 0 & shp.fort$exposure <= 0.01] <- 0
shp.fort$exposure_cat[shp.fort$exposure > 0.01 & shp.fort$exposure <= 0.02] <- 1
shp.fort$exposure_cat[shp.fort$exposure > 0.02 & shp.fort$exposure <= 0.03] <- 2
shp.fort$exposure_cat[shp.fort$exposure > 0.03 & shp.fort$exposure <= 0.04] <- 3
shp.fort$exposure_cat[shp.fort$exposure > 0.04 & shp.fort$exposure <= 0.05] <- 4
shp.fort$exposure_cat[shp.fort$exposure > 0.05 & shp.fort$exposure <= 0.06] <- 5
shp.fort$exposure_cat <- factor(shp.fort$exposure_cat, levels=0:5, 
                                labels = c("0-1", ">1-2", ">2-3", ">3-4", ">4-5", ">5-6"))

ggplot() +
  geom_polygon(data=shp.fort, aes(x=long, y=lat, group=group, alpha=exposure_cat),
               col="gray20", fill="cornflowerblue") + 
  labs(alpha = "Exposure\n(share of AB60 licenses)") +
  theme(axis.title=element_blank(), 
        axis.text=element_blank(), 
        axis.ticks=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill = NA),
        legend.position="bottom",
        legend.title = element_text(size=14, face = "bold", hjust=1),
        legend.text = element_text(size=12))


# Figure 1b: Correlation with share of unauthorized populatoin
ggplot(data = d.county[d.county$year==2015,], 
       aes(x=undoc_share*100, y=exposure*100, size=population)) + 
  geom_smooth(method="lm", col="firebrick3", alpha=.2) + 
  geom_point(col="gray15", alpha=.7) + 
  geom_text_repel(aes(label=county),size=4) +  
  xlab("Share of unauthorized immigrants\n(percent)") + 
  ylab("Share of AB60 licenses\n(percent)") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))






### Figure 2: Insurance uptake over time
# Figure 2a: Share of vehicles without insurance (%)
agg <- aggregate(d.county$uninsured,
                 by = list(d.county$exposure_tercile, d.county$year),
                 FUN = mean, na.rm=T)
colnames(agg) <- c("exposure", "year", "uninsured")
agg$year <- as.numeric(as.character(agg$year))
agg$exposure <- factor(agg$exposure, levels=1:3, labels = c("Low", "Medium", "High"))

ggplot(data = agg[agg$year >= 2006 & agg$year <= 2018,], aes(x=year, y=uninsured*100, col=as.factor(exposure), lty=as.factor(exposure))) + 
  annotate("rect", xmin = 2014.5, xmax = 2018.5, ymin = 10, ymax = 22, alpha = .3, fill="cornflowerblue") +
  geom_line(size=1.3) + 
  coord_cartesian(xlim=c(2006,2018.5), ylim=c(10, 22)) +  
  scale_x_continuous(limits=c(2006,2018.5), breaks=seq(2006,2018,1)) +
  scale_color_manual(values=c("gray60", "gray55", "gray10"), name="Exposure") +
  scale_linetype_manual(values=c(3,2,1), name="Exposure") +
  xlab("Year") + 
  ylab("Share of uninsured vehicles\n(percent)") +
  theme_classic() +
  theme(legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))




# Figure 2b: Uninsured motorist BI claims (per 1,000 vehicles)
agg <- aggregate(d.county$ClaimsNormStandUMBIBasic,
                 by = list(d.county$exposure_tercile, d.county$year),
                 FUN = mean, na.rm=T)
colnames(agg) <- c("exposure", "year", "ClaimsNormStandUMBIBasic")
agg$year <- as.numeric(as.character(agg$year))
agg$exposure <- factor(agg$exposure, levels=1:3, labels = c("Low", "Medium", "High"))

ggplot(data = agg[agg$year <= 2018,], 
       aes(x=year, y=ClaimsNormStandUMBIBasic, col=as.factor(exposure), lty=as.factor(exposure))) + 
  annotate("rect", xmin = 2014.5, xmax = 2018.5, ymin = 0, ymax = 5, alpha = .3, fill="cornflowerblue") +
  geom_line(size=1.3) + 
  coord_cartesian(xlim=c(2006,2018.5), ylim=c(0, 5)) +  
  scale_x_continuous(limits=c(2006,2018.5), breaks=seq(2006,2018,1)) +
  scale_color_manual(values=c("gray60", "gray55", "gray10"), name="Exposure") +
  scale_linetype_manual(values=c(3,2,1), name="Exposure") +
  xlab("Year") + 
  ylab("Claims for uninsured motorist BI\n(per 1,000 Vehicles)") +
  theme_classic() +
  theme(legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))



### Table 2: AB60 did not affect insurance uptake
m1 <- felm(uninsured ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(uninsured ~ exposure*treat | county + year | 0 | county, data=d.county)
m3 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(ClaimsNormStandUMBIBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Share ininsured", "UMBI claims"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES")))



### Figure 3: Insurance premiums over time
# Figure 3a: Average premium for bodily injury
agg <- aggregate(d.county$AvgPremiumStandBIBasic,
                 by = list(d.county$exposure_tercile, d.county$year),
                 FUN = mean, na.rm=T)
colnames(agg) <- c("exposure", "year", "AvgPremiumStandBIBasic")
agg$year <- as.numeric(as.character(agg$year))
agg$exposure <- factor(agg$exposure, levels=1:3, labels = c("Low", "Medium", "High"))

ggplot(data = agg[agg$year <= 2018,], 
       aes(x=year, y=AvgPremiumStandBIBasic, col=as.factor(exposure), 
           lty=as.factor(exposure))) + 
  annotate("rect", xmin = 2014.5, xmax = 2018.5, ymin = 130, ymax = 200, alpha = .3, fill="cornflowerblue") +
  geom_line(size=1.3) + 
  coord_cartesian(xlim=c(2006,2018.5), ylim=c(130,200)) +  
  scale_x_continuous(limits=c(2006,2018.5), breaks=seq(2006,2018,1)) +
  scale_color_manual(values=c("gray60", "gray55", "gray10"), name="Exposure") +
  scale_linetype_manual(values=c(3,2,1), name="Exposure") +
  xlab("Year") + 
  ylab("Average monthly premiums for BI coverage\n(US$)") +
  theme_classic() +
  theme(legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))


# Figure 3b: Average premium for inexperienced drivers
agg <- aggregate(d.county$premium1both,
                 by = list(d.county$exposure_tercile, d.county$year),
                 FUN = mean, na.rm=T)
colnames(agg) <- c("exposure", "year", "premium1both")
agg$year <- as.numeric(as.character(agg$year))
agg$exposure <- factor(agg$exposure, levels=1:3, labels = c("Low", "Medium", "High"))

ggplot(data = agg[!is.na(agg$premium1both),], 
       aes(x=year, y=premium1both, col=as.factor(exposure), lty=as.factor(exposure))) + 
  annotate("rect", xmin = 2014.5, xmax = 2018.5, ymin = 400, ymax = 1150, alpha = .3, fill="cornflowerblue") +
  geom_line(size=1.3) + 
  coord_cartesian(xlim=c(2006,2018.5), ylim=c(400, 1150)) +  
  scale_x_continuous(limits=c(2006,2018.5), breaks=seq(2006,2018,1)) +
  scale_color_manual(values=c("gray60", "gray55", "gray10"), name="Exposure") +
  scale_linetype_manual(values=c(3,2,1), name="Exposure") +
  xlab("Year") + 
  ylab("Average annual premiums\n(US$)") +
  theme_classic() +
  theme(legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))


# Figure 3c: Average premium for experienced drivers
agg <- aggregate(d.county$premium7both,
                 by = list(d.county$exposure_tercile, d.county$year),
                 FUN = mean, na.rm=T)
colnames(agg) <- c("exposure", "year", "premium7both")
agg$year <- as.numeric(as.character(agg$year))
agg$exposure <- factor(agg$exposure, levels=1:3, labels = c("Low", "Medium", "High"))

ggplot(data = agg[!is.na(agg$premium7both),], aes(x=year, y=premium7both, col=as.factor(exposure), lty=as.factor(exposure))) + 
  annotate("rect", xmin = 2014.5, xmax = 2018.5, ymin = 400, ymax = 1150, alpha = .3, fill="cornflowerblue") +
  geom_line(size=1.3) + 
  coord_cartesian(xlim=c(2006,2018.5), ylim=c(400, 1150)) +  
  scale_x_continuous(limits=c(2006,2018.5), breaks=seq(2006,2018,1)) +
  scale_color_manual(values=c("gray60", "gray55", "gray10"), name="Exposure") +
  scale_linetype_manual(values=c(3,2,1), name="Exposure") +
  xlab("Year") + 
  ylab("Average annual premiums\n(US$)") +
  theme_classic() +
  theme(legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))





### Table 3: AB60 did not affect insurance premiums
m1 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(AvgPremiumStandBIBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county)
m3 <- felm(premium1both_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(premium1both_log ~ exposure*treat | county + year | 0 | county, data=d.county)
m5 <- felm(premium7both_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m6 <- felm(premium7both_log ~ exposure*treat | county + year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4, m5, m6,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Bodily injury", "Inexperienced drivers", "Experienced drivers"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))






### Table 4: AB60 did not affect insurance premiums 2: Triple differences
m1 <- felm(premium_log ~ as.factor(exposure_tercile)*treat*treat2 | county + year | 0 | county, data=d.aps)
m2 <- felm(premium_log ~ exposure*treat*treat2  | county + year | 0 | county, data=d.aps)

stargazer(m1, m2,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Insurance premium"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Inexperienced",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Medium exposure x Inexperienced",
                             "High exposure x Inexperienced",                             
                             "Post-AB60 x Exposure (continuous)",
                             "Exposure (continuous) x Inexperienced", 
                             "Post-AB60 x Inexperienced",
                             "Post-AB60 x Medium exposure x Inexperienced",
                             "Post-AB60 x High exposure x Inexperienced",
                             "Post-AB60 x Exposure (continuous) x Inexperienced"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES"),
            c("Year-FEs?", "YES", "YES")))



### Figure 4: Power analysis
# note: see separate STATA do-file on how the power estimates were derived.

# Figure 4a: Share uninsured vehicles
ggplot() + 
  geom_hline(yintercept = 80, lty=2, size=.5) + 
  geom_point(data = power_uptake_uninsured, aes(x = effect_size, y = significant*100), size=1.5) +
  geom_line(data = power_uptake_uninsured, aes(x = effect_size, y = significant*100), size=1) +
  scale_y_continuous(limits=c(0,100), breaks = seq(0,100,20)) + 
  scale_color_manual(values=c("gray75", "gray25"), name="Exposure") + 
  xlab("Pseudo-Treatment Effect") + 
  ylab("Power") + 
  theme_classic() + 
  theme(plot.title = element_text(size=16, face="bold", hjust=.5),
        axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=12),
        legend.position="bottom")

# Figure 4b: UMBI claims
ggplot() + 
  geom_hline(yintercept = 80, lty=2, size=.5) + 
  geom_point(data = power_uptake_umbi, aes(x = effect_size, y = significant*100), size=1.5) +
  geom_line(data = power_uptake_umbi, aes(x = effect_size, y = significant*100), size=1) +
  scale_y_continuous(limits=c(0,100), breaks = seq(0,100,20)) + 
  scale_color_manual(values=c("gray75", "gray25"), name="Exposure") + 
  xlab("Pseudo-Treatment Effect") + 
  ylab("Power") + 
  theme_classic() + 
  theme(plot.title = element_text(size=16, face="bold", hjust=.5),
        axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=12),
        legend.position="bottom")

# Figure 4c: Bodily injury
ggplot() + 
  geom_hline(yintercept = 80, lty=2, size=.5) + 
  geom_point(data = power_premiums_bi, aes(x = effect_size, y = significant*100), size=1.5) +
  geom_line(data = power_premiums_bi, aes(x = effect_size, y = significant*100), size=1) +
  scale_y_continuous(limits=c(0,100), breaks = seq(0,100,20)) + 
  scale_color_manual(values=c("gray75", "gray25"), name="Exposure") + 
  xlab("Pseudo-Treatment Effect") + 
  ylab("Power") + 
  theme_classic() + 
  theme(plot.title = element_text(size=16, face="bold", hjust=.5),
        axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=12),
        legend.position="bottom")

# Figure 4d: Inexperienced drivers
ggplot() + 
  geom_hline(yintercept = 80, lty=2, size=.5) + 
  geom_point(data = power_premiums_inexperienced, aes(x = effect_size, y = significant*100), size=1.5) +
  geom_line(data = power_premiums_inexperienced, aes(x = effect_size, y = significant*100), size=1) +
  scale_y_continuous(limits=c(0,100), breaks = seq(0,100,20)) + 
  scale_color_manual(values=c("gray75", "gray25"), name="Exposure") + 
  xlab("Pseudo-Treatment Effect") + 
  ylab("Power") + 
  theme_classic() + 
  theme(plot.title = element_text(size=16, face="bold", hjust=.5),
        axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=12),
        legend.position="bottom")

# Figure 4e: Experienced drivers
ggplot() + 
  geom_hline(yintercept = 80, lty=2, size=.5) + 
  geom_point(data = power_premiums_experienced, aes(x = effect_size, y = significant*100), size=1.5) +
  geom_line(data = power_premiums_experienced, aes(x = effect_size, y = significant*100), size=1) +
  scale_y_continuous(limits=c(0,100), breaks = seq(0,100,20)) + 
  scale_color_manual(values=c("gray75", "gray25"), name="Exposure") + 
  xlab("Pseudo-Treatment Effect") + 
  ylab("Power") + 
  theme_classic() + 
  theme(plot.title = element_text(size=16, face="bold", hjust=.5),
        axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=12),
        legend.position="bottom")



### Table 5: The effect of AB60 on additional outcomes
m1 <- felm(autos_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(autos_log ~ exposure*treat | county + year | 0 | county, data=d.county)
m3 <- felm(PremiumStandBIBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(PremiumStandBIBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county)
m5 <- felm(ExposuresStandBIBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m6 <- felm(ExposuresStandBIBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county)
m7 <- felm(claimfrequency ~ as.factor(exposure_tercile)*treat | 
             county + year | 0 | county, data=d.county)
m8 <- felm(claimfrequency ~ exposure*treat | 
             county + year | 0 | county, data=d.county)

stargazer(m1,m2,m3,m4,m5,m6,m7,m8,
          df = F,
          omit.stat = c("rsq", "f"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Registered cars", "Premium volume", "Exposure volume", "Claim frequency"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES")))






# Appendix ----------------------------------------------------------------

######## A Summary Statistics ########

### Table A1: Summary statistics
# write function
summstats <- function(dv){
  # mean and sd for entire sample
  mean1 <- mean(d.county[,dv], na.rm=T)
  sd1 <- sd(d.county[,dv], na.rm=T)
  
  # mean and sd for pre-AB60 sample
  mean2 <- mean(d.county[as.numeric(as.character(d.county$year)) < 2015,dv], na.rm=T)
  sd2 <- sd(d.county[as.numeric(as.character(d.county$year)) < 2015,dv], na.rm=T)
  
  # mean and sd for post-AB60 sample
  mean3 <- mean(d.county[as.numeric(as.character(d.county$year)) >= 2015,dv], na.rm=T)
  sd3 <- sd(d.county[as.numeric(as.character(d.county$year)) >= 2015,dv], na.rm=T)  
  
  # create output
  out <- c(mean1, sd1, mean2, sd2, mean3, sd3)
  
  # return output
  return(out)
}

# table
stats <- data.frame(rbind(summstats(dv = "exposure"),summstats(dv = "uninsured"), 
                          summstats(dv = "ClaimsNormStandUMBIBasic"), summstats(dv = "ClaimsNormStandUMBIBasic_log"),
                          summstats(dv = "AvgPremiumStandBIBasic"), summstats(dv = "AvgPremiumStandBIBasic_log"), 
                          summstats(dv = "premium1both"), summstats(dv = "premium1both_log"), 
                          summstats(dv = "premium7both"), summstats(dv = "premium7both_log"),
                          summstats(d = "autos"), summstats(d = "autos_log"), 
                          summstats(d = "claimfrequency"),
                          summstats(d = "PremiumStandBIBasic"), summstats(d = "PremiumStandBIBasic_log"),
                          summstats(d = "ExposuresStandBIBasic"), summstats(d = "ExposuresStandBIBasic_log"),
                          summstats(d = "ClaimsNormNStandBIBasic"), summstats(d = "ClaimsNormNStandBIBasic_log"),
                          summstats(dv = "exposureNEW"),
                          summstats(dv = "undoc_share"),
                          summstats(dv = "UMBI_BI_ratio"), 
                          summstats(dv = "ClaimsNormStandUMPDBasic"), 
                          summstats(dv = "AvgPremiumStandPDBasic"), summstats(dv = "AvgPremiumStandPDBasic_log"),
                          summstats(dv = "unemprate"), summstats(dv = "popdens_log"),
                          summstats(dv = "income_log")))

colnames(stats) <- c("Mean", "SD", "Mean (Pre-AB60)", "SD (Pre-AB60)", "Mean (Post-AB60)", "SD (Post-AB60)")
rownames(stats) <- c("Share of AB60 licenses (2015)",
                     "Share of uninsured vehicles",
                     "Claims for uninsured motorist BI (per 1,000 vehicles)",
                     "Claims for uninsured motorist BI (per 1,000 vehicles; log)",
                     "Average monthly premium for bodily injury",
                     "Average monthly premium for bodily injury (log)",
                     "Annual premium for inexperienced drivers",
                     "Annual premium for inexperienced drivers (log)",
                     "Annual premium for experienced drivers",
                     "Annual premium for experienced drivers (log)",
                     "Registered cars", "Registered cars (log)",
                     "Claim frequency",
                     "Premium volume", "Premium volume (log)", 
                     "Exposure volume", "Exposure volume (log)", 
                     "Claims volume", "Claims volume (log)",
                     "Share of AB60 licenses (2015-16)",
                     "Share unauthorized immigrants",
                     "UMBI/BI ratio",
                     "Claims for uninsured motorist PD (per 1,000 vehicles)",
                     "Average monthly premium for property damage",
                     "Average monthly premium for property damage (log)",
                     "Unemployment rate", "Population density (log)", "Household income (log)")

print(xtable(stats, digits=3, align = "rcccccc", caption="Summary Statistics"), caption.placement = "top")






######## B Event Study Plots ########

### Figure A1: Event study plots for insurance uptake
# Figure A1a: Share of vehicles without insurance
m1 <- felm(uninsured ~ exposure*as.factor(yearEVENT) | county | 0 | county, data=d.county)

coefs <- data.frame(matrix(nrow=13, ncol=3))
coefs[,1] <- c(2006:2018)
coefs[c(1:8,10:13),2:3] <- coef(summary(m1))[14:25,1:2]
coefs[9,2:3] <- 0
colnames(coefs) <- c("year", "b", "se")

ggplot() + 
  annotate("rect", xmin=0.5, xmax=4.5, ymin=-Inf, ymax=Inf, fill="cornflowerblue", alpha=.4) + 
  geom_hline(yintercept = 0, lty=2) + 
  geom_point(data=coefs, aes(x=year-2014, y=b)) + 
  geom_errorbar(data=coefs, aes(x=year-2014, ymin=b-1.96*se, ymax=b+1.96*se), width=0) +
  scale_x_continuous(breaks=-8:4) + 
  xlab("Years before/after implementation") + 
  ylab("Estimated coefficient") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))


# Figure A1b: UMBI claims
m2 <- felm(ClaimsNormStandUMBIBasic_log ~ exposure*as.factor(yearEVENT) | county | 0 | county, data=d.county)

coefs <- data.frame(matrix(nrow=13, ncol=3))
coefs[,1] <- c(2006:2018)
coefs[c(1:8,10:13),2:3] <- coef(summary(m2))[14:25,1:2]
coefs[9,2:3] <- 0
colnames(coefs) <- c("year", "b", "se")

ggplot() + 
  annotate("rect", xmin=0.5, xmax=4.5, ymin=-Inf, ymax=Inf, fill="cornflowerblue", alpha=.4) +  
  geom_hline(yintercept = 0, lty=2) + 
  geom_point(data=coefs, aes(x=year-2014, y=b)) + 
  geom_errorbar(data=coefs, aes(x=year-2014, ymin=b-1.96*se, ymax=b+1.96*se), width=0) +
  scale_x_continuous(breaks=-8:4) + 
  xlab("Years before/after implementation") + 
  ylab("Estimated coefficient") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))



### Figure A2: Event study plots for insurance premiums
# Figure A2a: Average premium for bodily injury
m1 <- felm(AvgPremiumStandBIBasic_log ~ exposure*as.factor(yearEVENT) | county | 0 | county, data=d.county)

coefs <- data.frame(matrix(nrow=13, ncol=3))
coefs[,1] <- 2006:2018
coefs[c(1:8,10:13),2:3] <- coef(summary(m1))[14:25,1:2]
coefs[9,2:3] <- 0
colnames(coefs) <- c("year", "b", "se")

ggplot() + 
  annotate("rect", xmin=0.5, xmax=4.5, ymin=-Inf, ymax=Inf, fill="cornflowerblue", alpha=.4) +  
  geom_hline(yintercept = 0, lty=2) + 
  geom_point(data=coefs, aes(x=year-2014, y=b)) + 
  geom_errorbar(data=coefs, aes(x=year-2014, ymin=b-1.96*se, ymax=b+1.96*se), width=0) +
  scale_x_continuous(breaks=-8:4) + 
  xlab("Years before/after implementation") + 
  ylab("Estimated coefficient") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))





# Figure A2b: Average premium for inexperienced drivers
m2 <- felm(premium1both_log ~ exposure*as.factor(yearEVENT) | county | 0 | county, data=d.county)

coefs <- data.frame(matrix(nrow=13, ncol=3))
coefs[,1] <- 2006:2018
coefs[c(1:3),2:3] <- coef(summary(m2))[13:15,1:2]
coefs[c(5:8,10:13),2:3] <- coef(summary(m2))[16:23,1:2]
coefs[9,2:3] <- 0
colnames(coefs) <- c("year", "b", "se")

ggplot() + 
  annotate("rect", xmin=0.5, xmax=4.5, ymin=-Inf, ymax=Inf, fill="cornflowerblue", alpha=.4) +  
  geom_hline(yintercept = 0, lty=2) + 
  geom_point(data=coefs, aes(x=year-2014, y=b)) + 
  geom_errorbar(data=coefs, aes(x=year-2014, ymin=b-1.96*se, ymax=b+1.96*se), width=0) +
  scale_x_continuous(breaks=-8:4) + 
  xlab("Years before/after implementation") + 
  ylab("Estimated coefficient") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))



# Figure A2c: Average premium for experienced drivers
m3 <- felm(premium7both_log ~ exposure*as.factor(yearEVENT) | county | 0 | county, data=d.county)

coefs <- data.frame(matrix(nrow=13, ncol=3))
coefs[,1] <- 2006:2018
coefs[c(1:3),2:3] <- coef(summary(m3))[13:15,1:2]
coefs[c(5:8,10:13),2:3] <- coef(summary(m3))[16:23,1:2]
coefs[9,2:3] <- 0
colnames(coefs) <- c("year", "b", "se")

ggplot() + 
  annotate("rect", xmin=0.5, xmax=4.5, ymin=-Inf, ymax=Inf, fill="cornflowerblue", alpha=.4) +  
  geom_hline(yintercept = 0, lty=2) + 
  geom_point(data=coefs, aes(x=year-2014, y=b)) + 
  geom_errorbar(data=coefs, aes(x=year-2014, ymin=b-1.96*se, ymax=b+1.96*se), width=0) +
  scale_x_continuous(breaks=-8:4) + 
  xlab("Years before/after implementation") + 
  ylab("Estimated coefficient") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))



# Figure A2d: Triple difference
m1 <- felm(premium_log ~ exposure*as.factor(yearEVENT)*treat2 | county + year | 0 | county, data=d.aps)

coefs <- data.frame(matrix(nrow=13, ncol=3))
coefs[,1] <- 2006:2018
coefs[c(1:3, 5:8, 10:13),2:3] <- coef(summary(m1))[37:47,1:2]
coefs[9,2:3] <- 0
colnames(coefs) <- c("year", "b", "se")

ggplot() + 
  annotate("rect", xmin=0.5, xmax=4.5, ymin=-Inf, ymax=Inf, fill="cornflowerblue", alpha=.4) + 
  geom_hline(yintercept = 0, lty=2) + 
  geom_point(data=coefs, aes(x=year-2014, y=b)) + 
  geom_errorbar(data=coefs, aes(x=year-2014, ymin=b-1.96*se, ymax=b+1.96*se), width=0) +
  scale_x_continuous(breaks=-8:4) + 
  xlab("Years before/after implementation") + 
  ylab("Estimated coefficient") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))




######## C Robustness Checks ########

##### C.1 Alternative measures of insurance uptake and premiums #####

### Table A2: AB60 did not affect insurance uptake and premiums (alternative outcomes)
m1 <- felm(UMBI_BI_ratio ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(UMBI_BI_ratio ~ exposure*treat | county + year | 0 | county, data=d.county)
m3 <- felm(ClaimsNormStandUMPDBasic ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(ClaimsNormStandUMPDBasic ~ exposure*treat | county + year | 0 | county, data=d.county)
m5 <- felm(AvgPremiumStandPDBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county)
m6 <- felm(AvgPremiumStandPDBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4, m5, m6,
          df = F,
          omit.stat = c("rsq", "f"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("UMBI/BI ratio", "UMPD claims", "Premiums for PD"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))


### Table A3: AB60 did not affect the UMBI/BI ratio (different measures of exposure)
m1 <- felm(UMBI_BI_ratio ~ as.factor(exposureNEW_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(UMBI_BI_ratio ~ exposureNEW*treat | county + year | 0 | county, data=d.county)
m3 <- felm(UMBI_BI_ratio ~ as.factor(exposure_mean)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(UMBI_BI_ratio ~ as.factor(exposure_median)*treat | county + year | 0 | county, data=d.county)
m5 <- felm(UMBI_BI_ratio ~ as.factor(exposure_tercileDICH)*treat | county + year | 0 | county, data=d.county)
m6 <- felm(UMBI_BI_ratio ~ as.factor(undoc_share_tercile)*treat | county + year | 0 | county, data=d.county)
m7 <- felm(UMBI_BI_ratio ~ undoc_share*treat | county + year | 0 | county, data=d.county)

stargazer(m1,m2,m3,m4,m5,m6,m7,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("UMBI/BI ratio"),
          covariate.labels=c("Medium exposure (2015/16)",
                             "High exposure (2015/16",
                             "Exposure (2015/16; continuous)",
                             "Exposure (mean)",
                             "Exposure (median)",
                             "Exposure (low/medium vs. high)",
                             "Medium exposure (unauthorized)",
                             "High exposure (unauthorized)",
                             "Exposure (unauthorized; continous)",                             
                             "Post-AB60",
                             "Post-AB60 x Medium exposure (2015/16)",
                             "Post-AB60 x High exposure (2015/16)",
                             "Post-AB60 x Exposure (2015/16; continuous)",
                             "Post-AB60 x Exposure (mean)",
                             "Post-AB60 x Exposure (median)",
                             "Post-AB60 x Exposure (low/medium vs. high)",
                             "Post-AB60 x Medium exposure (unauthorized)",
                             "Post-AB60 x High exposure (unauthorized)",
                             "Post-AB60 x Exposure (unauthorized; continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES", "YES")))




##### C.2 Regressions weighted by population #####



### Table A4: AB60 did not affect insurance uptake (weighted regression)
m1 <- felm(uninsured ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))
m2 <- felm(uninsured ~ exposure*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))
m3 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))
m4 <- felm(ClaimsNormStandUMBIBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))

stargazer(m1,m2,m3,m4,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Share uninsured", "UMBI claims"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES")))




### Table A5: AB60 did not affect insurance premiums (weighted regression)
m1 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))
m2 <- felm(AvgPremiumStandBIBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))
m3 <- felm(premium1both_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))
m4 <- felm(premium1both_log ~ exposure*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))
m5 <- felm(premium7both_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))
m6 <- felm(premium7both_log ~ exposure*treat | county + year | 0 | county, data=d.county, 
           weights=log(d.county$pop.mean))

stargazer(m1,m2,m3,m4,m5,m6,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Bodily injury", "Inexperienced drivers", "Experienced drivers"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"),
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))




### Table A6: AB60 did not affect insurance premiums 2 (weighted regression)
m1 <- felm(premium_log ~ as.factor(exposure_tercile)*treat*treat2 | county + year | 0 | county, data=d.aps, 
           weights=log(d.aps$pop.mean))
m2 <- felm(premium_log ~ exposure*treat*treat2 | county + year | 0 | county, data=d.aps, 
           weights=log(d.aps$pop.mean))

stargazer(m1,m2,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Insurance premium"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Inexperienced",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Medium exposure x Inexperienced",
                             "High exposure x Inexperienced",                             
                             "Post-AB60 x Exposure (continuous)",
                             "Exposure (continuous) x Inexperienced", 
                             "Post-AB60 x Inexperienced",
                             "Post-AB60 x Medium exposure x Inexperienced",
                             "Post-AB60 x High exposure x Inexperienced",
                             "Post-AB60 x Exposure (continuous) x Inexperienced"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES"),
            c("Year-FEs?", "YES", "YES")))



##### C.3 Adding control variables #####

### Table A7: AB60 did not affect insurance uptake (with controls)
m1 <- felm(uninsured ~ as.factor(exposure_tercile)*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)
m2 <- felm(uninsured ~ exposure*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)
m3 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposure_tercile)*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)
m4 <- felm(ClaimsNormStandUMBIBasic_log ~ exposure*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)

stargazer(m1,m2,m3,m4,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Share uninsured", "UMBI claims"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Unemployment rate",
                             "Population density (log)",
                             "Income p.c. (log)",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES")))





### Table A8: AB60 did not affect insurance premiums (with controls)
m1 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposure_tercile)*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)
m2 <- felm(AvgPremiumStandBIBasic_log ~ exposure*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)
m3 <- felm(premium1both_log ~ as.factor(exposure_tercile)*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)
m4 <- felm(premium1both_log ~ exposure*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)
m5 <- felm(premium7both_log ~ as.factor(exposure_tercile)*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)
m6 <- felm(premium7both_log ~ exposure*treat + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.county)

stargazer(m1,m2,m3,m4,m5,m6,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Bodily injury", "Inexperienced drivers", "Experienced drivers"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Unemployment rate",
                             "Population density (log)",
                             "Income p.c. (log)",                             
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))




### Table A9: AB60 did not affect insurance premiums 2 (with controls)
m1 <- felm(premium_log ~ as.factor(exposure_tercile)*treat*treat2 + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.aps)
m2 <- felm(premium_log ~ exposure*treat*treat2 + 
             unemprate + popdens_log + income_log | county + year | 0 | county, data=d.aps)

stargazer(m1,m2,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Insurance premium"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Inexperienced",
                             "Unemployment rate",
                             "Population density (log)",
                             "Income p.c. (log)",                                
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Medium exposure x Inexperienced",
                             "High exposure x Inexperienced",                             
                             "Post-AB60 x Exposure (continuous)",
                             "Exposure (continuous) x Inexperienced", 
                             "Post-AB60 x Inexperienced",
                             "Post-AB60 x Medium exposure x Inexperienced",
                             "Post-AB60 x High exposure x Inexperienced",
                             "Post-AB60 x Exposure (continuous) x Inexperienced"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES"),
            c("Year-FEs?", "YES", "YES")))




##### C.4 Alternative definitions of exposure #####

### Table A10: Correlation between different measures of exposure to AB60
correlations <- data.frame(matrix(nrow=8, ncol=8))
correlations[,1] <- c("", "", "2015","2015", "2015/16", "2015/16", "unauthorized", "unauthorized")
correlations[,2] <- c("", "", "continuous", "terciles", "continuous", "terciles", "continuous", "terciles")
correlations[1,] <- c("", "", "2015","2015", "2015/16", "2015/16", "unauthorized", "unauthorized")
correlations[2,] <- c("", "", "continuous", "terciles", "continuous", "terciles", "continuous", "terciles")

correlations[3,3] <- round(cor(d.county$exposure[d.county$year==2015], d.county$exposure[d.county$year==2015]), 2)
correlations[3,4] <- round(cor(d.county$exposure[d.county$year==2015], d.county$exposure_tercile[d.county$year==2015]), 2)
correlations[3,5] <- round(cor(d.county$exposure[d.county$year==2015], d.county$exposureNEW[d.county$year==2015]), 2)
correlations[3,6] <- round(cor(d.county$exposure[d.county$year==2015], d.county$exposureNEW_tercile[d.county$year==2015]), 2)
correlations[3,7] <- round(cor(d.county$exposure[d.county$year==2015], d.county$undoc_share[d.county$year==2015]), 2)
correlations[3,8] <- round(cor(d.county$exposure[d.county$year==2015], d.county$undoc_share_tercile[d.county$year==2015]), 2)
correlations[4,3] <- round(cor(d.county$exposure_tercile[d.county$year==2015], d.county$exposure[d.county$year==2015]), 2)
correlations[4,4] <- round(cor(d.county$exposure_tercile[d.county$year==2015], d.county$exposure_tercile[d.county$year==2015]), 2)
correlations[4,5] <- round(cor(d.county$exposure_tercile[d.county$year==2015], d.county$exposureNEW[d.county$year==2015]), 2)
correlations[4,6] <- round(cor(d.county$exposure_tercile[d.county$year==2015], d.county$exposureNEW_tercile[d.county$year==2015]), 2)
correlations[4,7] <- round(cor(d.county$exposure_tercile[d.county$year==2015], d.county$undoc_share[d.county$year==2015]), 2)
correlations[4,8] <- round(cor(d.county$exposure_tercile[d.county$year==2015], d.county$undoc_share_tercile[d.county$year==2015]), 2)
correlations[5,3] <- round(cor(d.county$exposureNEW[d.county$year==2015], d.county$exposure[d.county$year==2015]), 2)
correlations[5,4] <- round(cor(d.county$exposureNEW[d.county$year==2015], d.county$exposure_tercile[d.county$year==2015]), 2)
correlations[5,5] <- round(cor(d.county$exposureNEW[d.county$year==2015], d.county$exposureNEW[d.county$year==2015]), 2)
correlations[5,6] <- round(cor(d.county$exposureNEW[d.county$year==2015], d.county$exposureNEW_tercile[d.county$year==2015]), 2)
correlations[5,7] <- round(cor(d.county$exposureNEW[d.county$year==2015], d.county$undoc_share[d.county$year==2015]), 2)
correlations[5,8] <- round(cor(d.county$exposureNEW[d.county$year==2015], d.county$undoc_share_tercile[d.county$year==2015]), 2)
correlations[6,3] <- round(cor(d.county$exposureNEW_tercile[d.county$year==2015], d.county$exposure[d.county$year==2015]), 2)
correlations[6,4] <- round(cor(d.county$exposureNEW_tercile[d.county$year==2015], d.county$exposure_tercile[d.county$year==2015]), 2)
correlations[6,5] <- round(cor(d.county$exposureNEW_tercile[d.county$year==2015], d.county$exposureNEW[d.county$year==2015]), 2)
correlations[6,6] <- round(cor(d.county$exposureNEW_tercile[d.county$year==2015], d.county$exposureNEW_tercile[d.county$year==2015]), 2)
correlations[6,7] <- round(cor(d.county$exposureNEW_tercile[d.county$year==2015], d.county$undoc_share[d.county$year==2015]), 2)
correlations[6,8] <- round(cor(d.county$exposureNEW_tercile[d.county$year==2015], d.county$undoc_share_tercile[d.county$year==2015]), 2)
correlations[7,3] <- round(cor(d.county$undoc_share[d.county$year==2015], d.county$exposure[d.county$year==2015]), 2)
correlations[7,4] <- round(cor(d.county$undoc_share[d.county$year==2015], d.county$exposure_tercile[d.county$year==2015]), 2)
correlations[7,5] <- round(cor(d.county$undoc_share[d.county$year==2015], d.county$exposureNEW[d.county$year==2015]), 2)
correlations[7,6] <- round(cor(d.county$undoc_share[d.county$year==2015], d.county$exposureNEW_tercile[d.county$year==2015]), 2)
correlations[7,7] <- round(cor(d.county$undoc_share[d.county$year==2015], d.county$undoc_share[d.county$year==2015]), 2)
correlations[7,8] <- round(cor(d.county$undoc_share[d.county$year==2015], d.county$undoc_share_tercile[d.county$year==2015]), 2)
correlations[8,3] <- round(cor(d.county$undoc_share_tercile[d.county$year==2015], d.county$exposure[d.county$year==2015]), 2)
correlations[8,4] <- round(cor(d.county$undoc_share_tercile[d.county$year==2015], d.county$exposure_tercile[d.county$year==2015]), 2)
correlations[8,5] <- round(cor(d.county$undoc_share_tercile[d.county$year==2015], d.county$exposureNEW[d.county$year==2015]), 2)
correlations[8,6] <- round(cor(d.county$undoc_share_tercile[d.county$year==2015], d.county$exposureNEW_tercile[d.county$year==2015]), 2)
correlations[8,7] <- round(cor(d.county$undoc_share_tercile[d.county$year==2015], d.county$undoc_share[d.county$year==2015]), 2)
correlations[8,8] <- round(cor(d.county$undoc_share_tercile[d.county$year==2015], d.county$undoc_share_tercile[d.county$year==2015]), 2)

print(xtable(correlations, digits=2, align = "rrccccccc", caption="Correlation of different measures of exposure to AB60"), 
      caption.placement = "top", include.rownames=FALSE, include.colnames = FALSE)



### Figure A3: Relation between different estimates of exposure to AB60
ggplot(data=d.county[d.county$year==2015,], aes(x=exposure*100, y=exposureNEW*100, size=log(population))) + 
  geom_smooth(method="lm", col="firebrick3", alpha=.2) + 
  geom_point(col="gray15", alpha=.7) + 
  geom_text_repel(aes(label=county),size=4)  +
  xlab("Share of AB60 licenses (2015)") + 
  ylab("Share of AB60 licenses (2015-16)") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))



### Table A11: AB60 did not affect insurance uptake (exposure in 2015/16)
m1 <- felm(uninsured ~ as.factor(exposureNEW_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(uninsured ~ exposureNEW*treat | county + year | 0 | county, data=d.county)
m3 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposureNEW_tercile)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(ClaimsNormStandUMBIBasic_log ~ exposureNEW*treat | county + year | 0 | county, data=d.county)

stargazer(m1,m2,m3,m4,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Share uninsured", "UMBI claims"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES")))





### Table A12: AB60 did not affect insurance premiums (exposure in 2015/16)
m1 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposureNEW_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(AvgPremiumStandBIBasic_log ~ exposureNEW*treat | county + year | 0 | county, data=d.county)
m3 <- felm(premium1both_log ~ as.factor(exposureNEW_tercile)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(premium1both_log ~ exposureNEW*treat | county + year | 0 | county, data=d.county)
m5 <- felm(premium7both_log ~ as.factor(exposureNEW_tercile)*treat | county + year | 0 | county, data=d.county)
m6 <- felm(premium7both_log ~ exposureNEW*treat | county + year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4, m5, m6,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Bodily injury", "Inexperienced drivers", "Experienced drivers"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"),
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))




### Table A13: AB60 did not affect insurance premiums 2 (exposure in 2015/16)
m1 <- felm(premium_log ~ as.factor(exposureNEW_tercile)*treat*treat2 | county + year | 0 | county, data=d.aps)
m2 <- felm(premium_log ~ exposureNEW*treat*treat2 | county + year | 0 | county, data=d.aps)

stargazer(m1,m2,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Insurance premium"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Inexperienced",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Medium exposure x Inexperienced",
                             "High exposure x Inexperienced",                             
                             "Post-AB60 x Exposure (continuous)",
                             "Exposure (continuous) x Inexperienced", 
                             "Post-AB60 x Inexperienced",
                             "Post-AB60 x Medium exposure x Inexperienced",
                             "Post-AB60 x High exposure x Inexperienced",
                             "Post-AB60 x Exposure (continuous) x Inexperienced"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES"),
            c("Year-FEs?", "YES", "YES")))




### Figure A4: Relation between estimated exposure to AB60 and the share of the unauthorized population
ggplot(data = d.county[d.county$year==2015,], 
       aes(x=undoc_share*100, y=exposure*100, size=population)) + 
  geom_smooth(method="lm", col="firebrick3", alpha=.2) + 
  geom_point(col="gray15", alpha=.7) + 
  geom_text_repel(aes(label=county),size=4) +  
  xlab("Share of unauthorized immigrants\n(percent)") + 
  ylab("Share of AB60 licenses\n(percent)") + 
  theme_classic() + 
  theme(legend.position = "none",
        axis.title = element_text(size=16, face="bold"),
        axis.text = element_text(size=14),
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=13))









### Table A14: AB60 did not affect insurance uptake (share unauthorized)
m1 <- felm(uninsured ~ as.factor(undoc_share_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(uninsured ~ undoc_share*treat | county + year | 0 | county, data=d.county)
m3 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(undoc_share_tercile)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(ClaimsNormStandUMBIBasic_log ~ undoc_share*treat | county + year | 0 | county, data=d.county)

stargazer(m1,m2,m3,m4,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Share uninsured", "UMBI claims"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"),
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES")))





### Table A15: AB60 did not affect insurance premiums (share unauthorized)
m1 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(undoc_share_tercile)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(AvgPremiumStandBIBasic_log ~ undoc_share*treat | county + year | 0 | county, data=d.county)
m3 <- felm(premium1both_log ~ as.factor(undoc_share_tercile)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(premium1both_log ~ undoc_share*treat | county + year | 0 | county, data=d.county)
m5 <- felm(premium7both_log ~ as.factor(undoc_share_tercile)*treat | county + year | 0 | county, data=d.county)
m6 <- felm(premium7both_log ~ undoc_share*treat | county + year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4, m5, m6,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Bodily injury", "Inexperienced drivers", "Experienced drivers"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))




### Table A16: AB60 did not affect insurance premiums 2 (share unauthorized)
m1 <- felm(premium_log ~ as.factor(undoc_share_tercile)*treat*treat2 | county + year | 0 | county, data=d.aps)
m2 <- felm(premium_log ~ undoc_share*treat*treat2 | county + year | 0 | county, data=d.aps)

stargazer(m1,m2,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Insurance premium"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Inexperienced",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Medium exposure x Inexperienced",
                             "High exposure x Inexperienced",                             
                             "Post-AB60 x Exposure (continuous)",
                             "Exposure (continuous) x Inexperienced", 
                             "Post-AB60 x Inexperienced",
                             "Post-AB60 x Medium exposure x Inexperienced",
                             "Post-AB60 x High exposure x Inexperienced",
                             "Post-AB60 x Exposure (continuous) x Inexperienced"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES"),
            c("Year-FEs?", "YES", "YES")))




##### C.5 Different cutoffs to identify high-exposure counties #####

### Table A17: AB60 did not affect insurance uptake (different cutoffs)
m1 <- felm(uninsured ~ as.factor(exposure_mean)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(uninsured ~ as.factor(exposure_median)*treat | county + year | 0 | county, data=d.county)
m3 <- felm(uninsured ~ as.factor(exposure_tercileDICH)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposure_mean)*treat | county + year | 0 | county, data=d.county)
m5 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposure_median)*treat | county + year | 0 | county, data=d.county)
m6 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposure_tercileDICH)*treat | county + year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4, m5, m6,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.caption = "",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = c("Share uninsured", "UMBI claims"),
          covariate.labels=c("Exposure (mean)",
                             "Exposure (median)",
                             "Exposure (low/medium vs. high)",
                             "Post-AB60",
                             "Post-AB60 x Exposure (mean)",
                             "Post-AB60 x Exposure (median)",
                             "Post-AB60 x Exposure (low/medium vs. high)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))



### Table A18: AB60 did not affect insurance premiums (different cutoffs)
m1 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposure_mean)*treat | county + year | 0 | county, data=d.county)
m2 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposure_median)*treat | county + year | 0 | county, data=d.county)
m3 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposure_tercileDICH)*treat | county + year | 0 | county, data=d.county)
m4 <- felm(premium1both_log ~ as.factor(exposure_mean)*treat | county + year | 0 | county, data=d.county)
m5 <- felm(premium1both_log ~ as.factor(exposure_median)*treat | county + year | 0 | county, data=d.county)
m6 <- felm(premium1both_log ~ as.factor(exposure_tercileDICH)*treat | county + year | 0 | county, data=d.county)
m7 <- felm(premium7both_log ~ as.factor(exposure_mean)*treat | county + year | 0 | county, data=d.county)
m8 <- felm(premium7both_log ~ as.factor(exposure_median)*treat | county + year | 0 | county, data=d.county)
m9 <- felm(premium7both_log ~ as.factor(exposure_tercileDICH)*treat | county + year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4, m5, m6, m7, m8, m9,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Bodily injury", "Inexperienced drivers", "Experienced drivers"),
          covariate.labels=c("Exposure (mean)",
                             "Exposure (median)",
                             "Exposure (low/medium vs. high)",
                             "Post-AB60",
                             "Post-AB60 x Exposure (mean)",
                             "Post-AB60 x Exposure (median)",
                             "Post-AB60 x Exposure (low/medium vs. high)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES")))


### Table A19: AB60 did not affect insurance premiums 2 (different cutoffs)
m1 <- felm(premium_log ~ as.factor(exposure_mean)*treat*treat2 | county + year | 0 | county, data=d.aps)
m2 <- felm(premium_log ~ as.factor(exposure_median)*treat*treat2 | county + year | 0 | county, data=d.aps)
m3 <- felm(premium_log ~ as.factor(exposure_tercileDICH)*treat*treat2 | county + year | 0 | county, data=d.aps)

stargazer(m1, m2, m3,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Insurance premium"),
          covariate.labels=c("Exposure (mean)",
                             "Exposure (median)",
                             "Exposure (low/medium vs. high)",
                             "Post-AB60",
                             "Inexperienced",
                             "Post-AB60 x Exposure (mean)",
                             "Post-AB60 x Exposure (median)",
                             "Post-AB60 x Exposure (low/medium vs. high)",
                             "Exposure (mean) x Inexperienced",
                             "Exposure (median) x Inexperienced",  
                             "Exposure (low/medium vs. high) x Inexperienced",
                             "Post-AB60 x Inexperienced",
                             "Post-AB60 x Exposure (mean) x Inexperienced",
                             "Post-AB60 x Exposure (median) x Inexperienced",
                             "Post-AB60 x Exposure (low/medium vs. high) x Inexperienced"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES")))





##### C.6 Shorter pre-treatment period #####


### Table A20: AB60 did not affect insurance uptake (different period)
m1 <- felm(uninsured ~ exposure*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)
m2 <- felm(uninsured ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)
m3 <- felm(ClaimsNormStandUMBIBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)
m4 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)

stargazer(m1, m2, m3, m4,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Share uninsured", "UMBI claims"),
          covariate.labels=c("Exposure (continuous)",
                             "Medium exposure",
                             "High exposure",
                             "Post-AB60",
                             "Post-AB60 x Exposure (continuous)",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES")))



### Table A21: AB60 did not affect insurance premiums (different period)
m1 <- felm(AvgPremiumStandBIBasic_log ~ exposure*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)
m2 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)
m3 <- felm(premium1both_log ~ exposure*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)
m4 <- felm(premium1both_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)
m5 <- felm(premium7both_log ~ exposure*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)
m6 <- felm(premium7both_log ~ as.factor(exposure_tercile)*treat | county + year | 0 | county, data=d.county, subset=year >= 2012)

stargazer(m1, m2, m3, m4, m5, m6, 
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Bodily injury", "Inexperienced drivers", "Experienced drivers"),
          covariate.labels=c("Exposure (continuous)",
                             "Medium exposure",
                             "High exposure",
                             "Post-AB60",
                             "Post-AB60 x Exposure (continuous)",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))


### Table A22: AB60 did not affect insurance premiums 2 (different period)
m1 <- felm(premium_log ~ exposure*treat*treat2 | county + year | 0 | county, data=d.aps, subset=year >= 2012)
m2 <- felm(premium_log ~ as.factor(exposure_tercile)*treat*treat2 | county + year | 0 | county, data=d.aps, subset=year >= 2012)

stargazer(m1, m2,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Insurance premium"),
          covariate.labels=c("Exposure (continuous)",
                             "Medium exposure",
                             "High exposure",
                             "Post-AB60",
                             "Inexperienced",
                             "Post-AB60 x Exposure (continuous)",                             
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Exposure (continuous) x Inexperienced",                             
                             "Medium exposure x Inexperienced",
                             "High exposure x Inexperienced", 
                             "Post-AB60 x Inexperienced",
                             "Post-AB60 x Exposure (continuous) x Inexperienced",
                             "Post-AB60 x Medium exposure x Inexperienced",
                             "Post-AB60 x High exposure x Inexperienced"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))



##### C.7 Dropping the county-fixed effects #####

### Table A23: AB60 did not affect insurance uptake (no county-FEs)
m1 <- felm(uninsured ~ as.factor(exposure_tercile)*treat | year | 0 | county, data=d.county)
m2 <- felm(uninsured ~ exposure*treat | year | 0 | county, data=d.county)
m3 <- felm(ClaimsNormStandUMBIBasic_log ~ as.factor(exposure_tercile)*treat | year | 0 | county, data=d.county)
m4 <- felm(ClaimsNormStandUMBIBasic_log ~ exposure*treat | year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Share uninsured", "UMBI claims"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES")))





### Table A24: AB60 did not affect insurance premiums (no county-FEs)
m1 <- felm(AvgPremiumStandBIBasic_log ~ as.factor(exposure_tercile)*treat | year | 0 | county, data=d.county)
m2 <- felm(AvgPremiumStandBIBasic_log ~ exposure*treat | year | 0 | county, data=d.county)
m3 <- felm(premium1both_log ~ as.factor(exposure_tercile)*treat | year | 0 | county, data=d.county)
m4 <- felm(premium1both_log ~ exposure*treat | year | 0 | county, data=d.county)
m5 <- felm(premium7both_log ~ as.factor(exposure_tercile)*treat | year | 0 | county, data=d.county)
m6 <- felm(premium7both_log ~ exposure*treat | year | 0 | county, data=d.county)

stargazer(m1, m2, m3, m4, m5, m6,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Bodily injury", "Inexperienced crivers", "Experienced crivers"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Post-AB60 x Exposure (continuous)"),
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FEs?", "YES", "YES", "YES", "YES", "YES", "YES")))




### Table A25: AB60 did not affect insurance premiums 2 (no county-FEs)
m1 <- felm(premium_log ~ as.factor(exposure_tercile)*treat*treat2 | year | 0 | county, data=d.aps)
m2 <- felm(premium_log ~ exposure*treat*treat2  | year | 0 | county, data=d.aps)

stargazer(m1, m2,
          df = F,
          omit.stat = c("rsq", "f", "ser"),
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels.include = T,
          dep.var.labels  = c("Insurance premium"),
          covariate.labels=c("Medium exposure",
                             "High exposure",
                             "Exposure (continuous)",
                             "Post-AB60",
                             "Inexperienced",
                             "Post-AB60 x Medium exposure",
                             "Post-AB60 x High exposure",
                             "Medium exposure x Inexperienced",
                             "High exposure x Inexperienced",                             
                             "Post-AB60 x Exposure (continuous)",
                             "Exposure (continuous) x Inexperienced", 
                             "Post-AB60 x Inexperienced",
                             "Post-AB60 x Medium exposure x Inexperienced",
                             "Post-AB60 x High exposure x Inexperienced",
                             "Post-AB60 x Exposure (continuous) x Inexperienced"), 
          digits = 3,
          add.lines = list(
            c("County-FEs?", "YES", "YES"),
            c("Year-FEs?", "YES", "YES")))

